This script is similar to the ek_irrad_model but uses relative Hsat for analysis The amount of time in Hsat was divided by the total day length, which is defined by the amount of time when umols photons m-2 s-1 were greater than or equal to 1.

#load the various libraries
library(lme4)
library(lmerTest)
library(effects)
library(car)
library(MuMIn)
library (dplyr)
library(emmeans)
library(DHARMa)
library(performance)
library(patchwork)
library(rstatix)
#for plots and tables
library(ggplot2)
library(ggpubr)
library(forcats)
library(RColorBrewer)
library(tidyverse)
library(sjPlot)
library(sjmisc)
library(mmtable2)
library(gt)
library(purrr)
library(stringr)
library(tidyr)

Load the dataset

ek_irrad_data <- read.csv("input_data/mean_supersaturation_by_period.csv")

Assign various variables as factors

ek_irrad_data$Run <- as.factor(ek_irrad_data$Run)
ek_irrad_data$Temperature <- as.factor(ek_irrad_data$Temp...C.)
ek_irrad_data$Treatment <- as.factor(as.character(ek_irrad_data$Treatment))
#ek_irrad_data$deltaNPQ <- as.factor(ek_irrad_data$deltaNPQ)

ULVA____________________________________________________________

Subset the data by species and make new column of the treatments for the plots

ulva <- subset(ek_irrad_data, Species == "ul" & Treatment != 2.5)
ulva$treatment_graph[ulva$Treatment == 0] <- "1) 35ppt/0.5umol"
ulva$treatment_graph[ulva$Treatment == 1] <- "2) 35ppt/14umol" 
ulva$treatment_graph[ulva$Treatment == 2] <- "3) 28ppt/27umol" 
ulva$treatment_graph[ulva$Treatment == 3] <- "5) 18ppt/53umol" 
ulva$treatment_graph[ulva$Treatment == 4] <- "6) 11ppt/80umol"
ulva$treatment_graph[ulva$Treatment == 2.5] <- "4) 28ppt/53umol"

Make a histogram for ulva

hist(ulva$supersat_rel, main = paste("Ulva lactuca"), col = "olivedrab3", labels = TRUE)

ulva %>% ggplot(aes(supersat_avg)) +
        geom_histogram(binwidth=5, fill = "#5BB300", color = "black", size = 0.25, alpha = 0.85) +
        theme_bw()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.

Run model without interaction between the treatments and temperature Take RLC.Order out of the random effects because causing problems of singularity. R2 is same with or without (+ (1 | RLC.Order))

rel_hsat_model_ulva <- lmer(formula = supersat_rel ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva)

Make residual plots of the data

hist(resid(rel_hsat_model_ulva))

plot(resid(rel_hsat_model_ulva) ~ fitted(rel_hsat_model_ulva))

qqnorm(resid(rel_hsat_model_ulva))
qqline(resid(rel_hsat_model_ulva))

Check the performance of the model, get R2, plot the effects of model and make table

performance::check_model(rel_hsat_model_ulva)
## Variable `Component` is not in your data frame :/

r.squaredGLMM(rel_hsat_model_ulva)
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
##            R2m      R2c
## [1,] 0.1917698 0.822725
summary(rel_hsat_model_ulva)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: supersat_rel ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +  
##     (1 | RLC.Order)
##    Data: ulva
## 
## REML criterion at convergence: 1443.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8397 -0.4267  0.0763  0.5141  2.6191 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Plant.ID  (Intercept) 14.3670  3.790   
##  Run       (Intercept) 39.1517  6.257   
##  RLC.Order (Intercept)  0.9662  0.983   
##  Residual              15.3082  3.913   
## Number of obs: 240, groups:  Plant.ID, 96; Run, 7; RLC.Order, 6
## 
## Fixed effects:
##               Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)    74.4406     4.5814   5.2312  16.248 1.12e-05 ***
## Treatment1     -5.8265     5.3559   4.9914  -1.088   0.3264    
## Treatment2     -5.2599     5.3559   4.9914  -0.982   0.3712    
## Treatment3    -10.2494     5.3559   4.9914  -1.914   0.1139    
## Treatment4    -11.1703     5.3559   4.9914  -2.086   0.0915 .  
## Temperature27   0.9231     1.4722   5.3816   0.627   0.5563    
## Temperature30   1.6739     1.3225  24.8618   1.266   0.2173    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Trtmn2 Trtmn3 Trtmn4 Tmpr27
## Treatment1  -0.823                                   
## Treatment2  -0.823  0.989                            
## Treatment3  -0.823  0.989  0.989                     
## Treatment4  -0.823  0.989  0.989  0.989              
## Temperatr27 -0.151 -0.001 -0.001 -0.001 -0.001       
## Temperatr30 -0.145  0.000  0.000  0.000  0.000  0.460
plot(allEffects(rel_hsat_model_ulva))

tab_model(rel_hsat_model_ulva, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)
  supersat_rel
Predictors Estimates std. Error CI Statistic p df
(Intercept) 74.44 4.58 65.41 – 83.47 16.25 <0.001 229.00
Treatment [1] -5.83 5.36 -16.38 – 4.73 -1.09 0.278 229.00
Treatment [2] -5.26 5.36 -15.81 – 5.29 -0.98 0.327 229.00
Treatment [3] -10.25 5.36 -20.80 – 0.30 -1.91 0.057 229.00
Treatment [4] -11.17 5.36 -21.72 – -0.62 -2.09 0.038 229.00
Temperature [27] 0.92 1.47 -1.98 – 3.82 0.63 0.531 229.00
Temperature [30] 1.67 1.32 -0.93 – 4.28 1.27 0.207 229.00
Random Effects
σ2 15.31
τ00 Plant.ID 14.37
τ00 Run 39.15
τ00 RLC.Order 0.97
ICC 0.78
N Run 7
N Plant.ID 96
N RLC.Order 6
Observations 240
Marginal R2 / Conditional R2 0.192 / 0.823

Construct null model to perform likelihood ratio test REML must be FALSE

ulva_relhsat_treatment_null <- lmer(formula = supersat_rel ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
ulva_relhsat_model2 <- lmer(formula = supersat_rel ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
anova(ulva_relhsat_treatment_null, ulva_relhsat_model2)
## Data: ulva
## Models:
## ulva_relhsat_treatment_null: supersat_rel ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## ulva_relhsat_model2: supersat_rel ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                             npar    AIC    BIC  logLik deviance  Chisq Df
## ulva_relhsat_treatment_null    7 1544.8 1569.1 -765.39   1530.8          
## ulva_relhsat_model2           11 1481.6 1519.9 -729.81   1459.6 71.146  4
##                             Pr(>Chisq)    
## ulva_relhsat_treatment_null               
## ulva_relhsat_model2            1.3e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ulva_relhsat_temperature_null <- lmer(formula = supersat_rel ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
ulva_relhsat_model3 <- lmer(formula = supersat_rel ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
anova(ulva_relhsat_temperature_null, ulva_relhsat_model3)
## Data: ulva
## Models:
## ulva_relhsat_temperature_null: supersat_rel ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## ulva_relhsat_model3: supersat_rel ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                               npar    AIC    BIC  logLik deviance  Chisq Df
## ulva_relhsat_temperature_null    9 1479.0 1510.4 -730.53   1461.0          
## ulva_relhsat_model3             11 1481.6 1519.9 -729.81   1459.6 1.4286  2
##                               Pr(>Chisq)
## ulva_relhsat_temperature_null           
## ulva_relhsat_model3               0.4895

Plots

ulva %>% ggplot(aes(treatment_graph, supersat_rel)) + 
        geom_boxplot(size=0.5) + 
        geom_point(alpha = 0.75, size = 3, aes(color = Temperature), show.legend = TRUE) + 
        labs(x="Treatment (salinty/nitrate)", y= "Relative Hsat (%)", title= "C", subtitle = "Ulva lactuca") + 
        scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "18ppt/53umolN", "11ppt/80umolN")) + 
        ylim(20, 100) + stat_mean() + 
        geom_hline(yintercept=60, color = "red", size = 0.5, alpha = 0.5) +
        scale_color_manual(values = c("#295102", "#7CB950", "#BDE269")) +
        theme_bw() +
        theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
              plot.subtitle = element_text(face = "italic", size = 14, vjust = -18, hjust = 0.05))

Summarize the means for relative Hsat

ulva %>% group_by(Treatment) %>% summarise_at(vars(supersat_rel), list(mean = mean))
## # A tibble: 5 × 2
##   Treatment  mean
##   <fct>     <dbl>
## 1 0          75.3
## 2 1          70.1
## 3 2          70.7
## 4 3          65.7
## 5 4          64.7
ulva %>% group_by(Run) %>% summarise_at(vars(supersat_rel), list(mean = mean))
## # A tibble: 7 × 2
##   Run    mean
##   <fct> <dbl>
## 1 1      72.3
## 2 2      73.8
## 3 3      69.2
## 4 3.5    60.0
## 5 4      60.4
## 6 6      71.2
## 7 6.5    79.4
ulva %>% group_by(Run) %>% summarise_at(vars(day_length_avg), list(mean = mean))
## # A tibble: 7 × 2
##   Run    mean
##   <fct> <dbl>
## 1 1      654.
## 2 2      643.
## 3 3      626.
## 4 3.5    621.
## 5 4      611.
## 6 6      626.
## 7 6.5    633.

Add growth rate to this dataset

growth_rate <- read.csv("/Users/Angela/src/work/limu/algal_growth_photosynthesis/data_input/all_runs_growth_011723.csv")
growth_rate$Species <- as.factor(growth_rate$Species)
growth_rate$treatment <- as.factor(growth_rate$treatment)

Make a new column for weight change (difference final from initial)

growth_rate$growth_rate_percent <- (growth_rate$final.weight - growth_rate$Initial.weight) / growth_rate$Initial.weight * 100

Subset by species and add to current dataset

gr_ulva <- subset(growth_rate, Species == "Ul" & treatment != 2.5)
ulva$growth_rate <- round((gr_ulva$final.weight - gr_ulva$Initial.weight) / gr_ulva$Initial.weight * 100, digits = 2)

Plot a regression between the photosynthetic independent variables of interest and growth rate

ulva_growth_rel_hsat_graph <- ggplot(ulva, aes(x=supersat_rel, y=growth_rate)) + 
        geom_point(alpha = 0.5, size = 3, aes(color = Treatment), show.legend = TRUE) + 
        geom_smooth(method = "lm", col = "black") + theme_bw() + 
        labs(title = "A", subtitle = "Ulva lactuca", x = "Mean Rel. Hsat Time (%)", 
             y = "growth rate (%)") + stat_regline_equation(label.x = 75, label.y = 225) + 
        stat_cor(label.x = 75, label.y = 215) +
        theme(plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
              plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))
ulva_growth_rel_hsat_graph
## `geom_smooth()` using formula = 'y ~ x'

HYPNEA_____________________________________________________________________________________________________________ Subset for Hypnea

hypnea <- subset(ek_irrad_data, Species == "hm" & day1_rlc_time != "11:34:05")
hypnea$treatment_graph[hypnea$Treatment == 0] <- "1) 35ppt/0.5umol"
hypnea$treatment_graph[hypnea$Treatment == 1] <- "2) 35ppt/14umol" 
hypnea$treatment_graph[hypnea$Treatment == 2] <- "3) 28ppt/27umol" 
hypnea$treatment_graph[hypnea$Treatment == 3] <- "5) 18ppt/53umol" 
hypnea$treatment_graph[hypnea$Treatment == 4] <- "6) 11ppt/80umol"
hypnea$treatment_graph[hypnea$Treatment == 2.5] <- "4) 28ppt/53umol"

Make a histogram of the data for ulva

hist(hypnea$supersat_rel, main = paste("Hypnea musciformis"), col = "maroon", labels = TRUE)

hypnea %>% ggplot(aes(supersat_rel)) +
        geom_histogram(binwidth=5, fill = "#a8325e", color = "black", size = 0.25, alpha = 0.85) +
        theme_bw()

Run model

rel_hsat_model_hypnea  <- lmer(formula = supersat_rel ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea)

Plot residuals

hist(resid(rel_hsat_model_hypnea))

plot(resid(rel_hsat_model_hypnea) ~ fitted(rel_hsat_model_hypnea))

qqnorm(resid(rel_hsat_model_hypnea))
qqline(resid(rel_hsat_model_hypnea))

Check the performance of the model, get R2, make table and plot effects

performance::check_model(rel_hsat_model_hypnea)
## Variable `Component` is not in your data frame :/

r.squaredGLMM(rel_hsat_model_hypnea)
##            R2m       R2c
## [1,] 0.2345395 0.7173404
summary(rel_hsat_model_hypnea)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: supersat_rel ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +  
##     (1 | RLC.Order)
##    Data: hypnea
## 
## REML criterion at convergence: 1824.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3952 -0.4447  0.0417  0.5604  2.3878 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Plant.ID  (Intercept) 15.7258  3.9656  
##  Run       (Intercept) 24.7192  4.9718  
##  RLC.Order (Intercept)  0.8515  0.9228  
##  Residual              24.1774  4.9171  
## Number of obs: 286, groups:  Plant.ID, 144; Run, 8; RLC.Order, 6
## 
## Fixed effects:
##               Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)     65.012      3.732   5.286  17.419 7.14e-06 ***
## Treatment1      -7.659      4.361   5.043  -1.756 0.138852    
## Treatment2      -5.662      4.361   5.043  -1.298 0.250378    
## Treatment2.5   -10.804      6.224   4.561  -1.736 0.148734    
## Treatment3      -7.924      4.361   5.043  -1.817 0.128380    
## Treatment4     -10.498      4.363   5.056  -2.406 0.060613 .  
## Temperature27    4.541      1.410   4.683   3.221 0.025750 *  
## Temperature30    6.123      1.310  13.044   4.674 0.000432 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Trtmn2 Trt2.5 Trtmn3 Trtmn4 Tmpr27
## Treatment1  -0.811                                          
## Treatment2  -0.811  0.974                                   
## Treatmnt2.5 -0.568  0.486  0.486                            
## Treatment3  -0.811  0.974  0.974  0.486                     
## Treatment4  -0.810  0.973  0.973  0.486  0.973              
## Temperatr27 -0.179  0.001  0.001  0.000  0.001  0.001       
## Temperatr30 -0.174 -0.001 -0.001  0.000 -0.001  0.000  0.453
plot(allEffects(rel_hsat_model_hypnea))

tab_model(rel_hsat_model_hypnea, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)
  supersat_rel
Predictors Estimates std. Error CI Statistic p df
(Intercept) 65.01 3.73 57.66 – 72.36 17.42 <0.001 274.00
Treatment [1] -7.66 4.36 -16.24 – 0.93 -1.76 0.080 274.00
Treatment [2] -5.66 4.36 -14.25 – 2.92 -1.30 0.195 274.00
Treatment [2.5] -10.80 6.22 -23.06 – 1.45 -1.74 0.084 274.00
Treatment [3] -7.92 4.36 -16.51 – 0.66 -1.82 0.070 274.00
Treatment [4] -10.50 4.36 -19.09 – -1.91 -2.41 0.017 274.00
Temperature [27] 4.54 1.41 1.77 – 7.32 3.22 0.001 274.00
Temperature [30] 6.12 1.31 3.54 – 8.70 4.67 <0.001 274.00
Random Effects
σ2 24.18
τ00 Plant.ID 15.73
τ00 Run 24.72
τ00 RLC.Order 0.85
ICC 0.63
N Run 8
N Plant.ID 144
N RLC.Order 6
Observations 286
Marginal R2 / Conditional R2 0.235 / 0.717

Construct null model to perform likelihood ratio test REML must be FALSE

hypnea_relhsat_treatment_null <- lmer(formula = supersat_rel ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
hypnea_relhsat_model2 <- lmer(formula = supersat_rel ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
anova(hypnea_relhsat_treatment_null, hypnea_relhsat_model2)
## Data: hypnea
## Models:
## hypnea_relhsat_treatment_null: supersat_rel ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## hypnea_relhsat_model2: supersat_rel ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                               npar    AIC    BIC  logLik deviance  Chisq Df
## hypnea_relhsat_treatment_null    7 1887.1 1912.7 -936.53   1873.1          
## hypnea_relhsat_model2           12 1870.3 1914.2 -923.17   1846.3 26.723  5
##                               Pr(>Chisq)    
## hypnea_relhsat_treatment_null               
## hypnea_relhsat_model2          6.458e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hypnea_relhsat_temperature_null <- lmer(formula = supersat_rel ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
hypnea_relhsat_model3 <- lmer(formula = supersat_rel ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
anova(hypnea_relhsat_temperature_null, hypnea_relhsat_model3)
## Data: hypnea
## Models:
## hypnea_relhsat_temperature_null: supersat_rel ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## hypnea_relhsat_model3: supersat_rel ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                                 npar    AIC    BIC  logLik deviance  Chisq Df
## hypnea_relhsat_temperature_null   10 1886.6 1923.1 -933.28   1866.6          
## hypnea_relhsat_model3             12 1870.3 1914.2 -923.17   1846.3 20.213  2
##                                 Pr(>Chisq)    
## hypnea_relhsat_temperature_null               
## hypnea_relhsat_model3            4.081e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plots

hypnea %>% ggplot(aes(treatment_graph, supersat_rel)) + 
        geom_boxplot(size=0.5) + 
        geom_point(alpha = 0.75, size = 3, aes(color = Temperature), show.legend = TRUE) + 
        labs(x="Treatment (salinty/nitrate)", y= "Relative Hsat (%)", title= "D", subtitle = "Hypnea musciformis") + 
        scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "28ppt/53umolN", "18ppt/53umolN", "11ppt/80umolN")) + 
        ylim(20, 100) + stat_mean() + 
        geom_hline(yintercept=60, color = "red", size = 0.5, alpha = 0.5) +
        scale_color_manual(values = c("#9C0627", "#BB589F", "#F4B4E2")) +
        theme_bw() +
        theme(legend.position = c(0.88,0.88), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
              plot.subtitle = element_text(face = "italic", size = 14, vjust = -18, hjust = 0.05))

For regression with growth Hypnea remove hm6-4 on 11/12 that had no d9 RLC (final weight 0.1017) and hm6-4 on 10/29/21 because it was white and also looked dead

gr_hypnea <- subset(growth_rate, Species == "Hm" & final.weight != 0.1017 & growth_rate_percent > -87.96837)
hypnea$growth_rate <- round((gr_hypnea$final.weight - gr_hypnea$Initial.weight) / gr_hypnea$Initial.weight * 100, digits = 2)

Plot a regression between the photosynthetic independent variables of interest and growth rate

hypnea_growth_rel_hsat_graph <- ggplot(hypnea, aes(x=supersat_rel, y=growth_rate)) + 
        geom_point(alpha = 0.5, size = 3, aes(color = Treatment), show.legend = TRUE) + 
        geom_smooth(method = "lm", col = "black") + theme_bw() + 
        labs(title = "B", subtitle = "Hypnea musciformis", x = "Mean Rel. Hsat (%)", y = "growth rate (%)") + 
        stat_regline_equation(label.x = 65, label.y = 200) + 
        stat_cor(label.x = 65, label.y = 215) +
        theme(plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
              plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))
hypnea_growth_rel_hsat_graph
## `geom_smooth()` using formula = 'y ~ x'

Summarize the means for relative Hsat

hypnea %>% group_by(Treatment) %>% summarise_at(vars(supersat_rel), list(mean = mean))
## # A tibble: 6 × 2
##   Treatment  mean
##   <fct>     <dbl>
## 1 0          68.6
## 2 1          61.7
## 3 2          63.6
## 4 2.5        57.8
## 5 3          61.4
## 6 4          58.7
hypnea %>% group_by(Run) %>% summarise_at(vars(supersat_rel), list(mean = mean))
## # A tibble: 8 × 2
##   Run    mean
##   <fct> <dbl>
## 1 1      64.5
## 2 2      67.8
## 3 3      60.2
## 4 3.5    54.2
## 5 4      56.0
## 6 7      57.8
## 7 8      67.0
## 8 9      70.1
hypnea %>% group_by(Run) %>% summarise_at(vars(day_length_avg), list(mean = mean))
## # A tibble: 8 × 2
##   Run    mean
##   <fct> <dbl>
## 1 1      654.
## 2 2      645.
## 3 3      626.
## 4 3.5    621.
## 5 4      611.
## 6 7      688.
## 7 8      635.
## 8 9      636.